////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//Implementation for the Main method
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#include "StdAfx.h"


/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

void defaultOutput()
{
		printf("\n");
		printf("#########################################################\n");
		printf("## 411 Homework: Numerical Analysis                    ##\n");
		printf("## Brian Chen, Oct 20, 2012                            ##\n");
		printf("#########################################################\n");
		printf("\n");
		printf("Usage:\n");
		printf("\n");
		printf("invert \n");
		printf("With no arguments, invert returns this usage statement.\n");
		printf("\n");
		printf("invert <filename>\n");
		printf("Prints out two matrix inverses for the provided matrix.\n");
		printf(" -- the first inverse is computed with simple inversion.\n");
		printf(" -- the second is computed with partial pivoting\n");
		printf("\n");
		printf("invert -format\n");
		printf("Describes the file format for the input file.\n");
		printf("\n");
		printf("#########################################################\n");
		printf("\n");
		exit(0);
}

void fileFormat()
{
		printf("\n");
		printf("#########################################################\n");
		printf("## 411 Homework: Numerical Analysis                    ##\n");
		printf("## Brian Chen, Oct 20, 2012                            ##\n");
		printf("#########################################################\n");
		printf("\n");
		printf("File format for a 3rd order matrix (input):\n");
		printf(" - numbers below are space deliminated.  no tabs.\n");
		printf(" - There are three lines, each with 3 numbers.\n");
		printf(" - They will be parsed as doubles.\n");
		printf(" - Comment lines start with #.  Comment lines cannot.\n");
		printf("   interrupt the 3 lines with the matrix, but can come.\n");
		printf("   before or after.\n");
		printf("\n");
		printf("---- Matrix file below, not including this line ----\n");
		printf("# comment\n");
		printf("# comment\n");
		printf("1.2 3.3 4.4\n");
		printf("1.112 2.24 0\n");
		printf("3.4 2.2 98.23452\n");
		printf("#\n");
		printf("#\n");
		printf("---- Matrix file above, not including this line ----\n");
		printf("\n");
		printf("#########################################################\n");
		printf("\n");
		exit(0);
}


/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

//Inside the result, the result is to be stored as:
//i=0 i=1 i=2
//[0] [3] [6]  //j=0
//[1] [4] [7]  //j=1
//[2] [5] [8]  //j=2
//
//where the above 3x3 matrix is as you would write it on paper.
//
//
//
//


//This method swaps two double variables
void swap(double &a, double &b)
{
	double tmp;
	tmp = a;
	a = b;
	b = tmp;
}

//This method calculates the inverse of a matrix by using the gaussian elimination method
double ** simpleMatrixInversion( double ** matrix ){
	int i = 0;
	int j = 0;
	int k = 0;
	//insert your code here.//

	//allocate the matrix for calculation.
	double temp_matrix[3][3];

	//assign the matrix values that from the functional parameter matrix
	for(i = 0; i < 3; i ++)
	{
		for(j = 0; j < 3; j++)
		{
			temp_matrix[i][j] = matrix[i][j];
		}
	}

	
	//declare and allocate the result
	//initialize the result matrix as an identity matrix
	double ** result = new double*[3];
	for(i = 0; i<3; i++){
		result[i] = new double[3];
		for(j = 0; j<3; j++)
			result[i][j] = (i == j) ? 1.0 : 0.0;
	}

	//calculate the inverse of the matrix
	//here we do row transformation for the larger matrix of (temp_matrix, result)
	//when temp_matrix becomes identity matrix, the result matrix becomes the inverse of the original
	//temp_matrix
        
	for(i = 0; i<3; i++)
	{
		//if line i contains a zero in the diagonal place, i.e. place [i][i]
		//then we search each following lines to find an alternative line which does not contain
		//a zero in the place [alternative][i]
		int alternative = i;
		while(alternative < 3 && temp_matrix[alternative][i] == 0)
			++alternative;

		//row i and all rows below row  i contain zero in the column i, therefore, 
		//it can be proved that such a matrix do not have an inverse.
		if(alternative == 3)
		{
			printf("This matrix is singular and do not have an inverse!\n");
			exit(0);
		}

		//if the program runs here, it means that a non-zero alternative exists
		//if this alternative line different from line i
		//we swap each element in both rows
		if(alternative != i)
		{
			for(k = i; k<3; k++)
				swap(temp_matrix[alternative][k], temp_matrix[i][k]);
			for(k = 0; k<3; k++)
				swap(result[alternative][k], result[i][k]);
		}

		//for each row down from the ith row
                for(j = i + 1; j<3; j++)
                {
			if(temp_matrix[j][i] != 0)
			{
				//here we gurantee that matrix[i][i] is not zero
				double multiple = temp_matrix[j][i] / temp_matrix[i][i];
				for(k = i; k<3; k++)
                                	 temp_matrix[j][k] -= temp_matrix[i][k] * multiple;
				for(k = 0; k<3; k++)
					 result[j][k] -= result[i][k] * multiple;
			}
			else
			{
				//the target row has been zero, therefore, 
				//we just skip it and move on to the next row
				continue;
			}

		}	
	}


	//we do the substition process for each line from the largest line to the smallest. i.e., 0
	for(i = 2; i >= 0; i--)
	{
		//try to make the ith diagonal place 1
		double multiple = temp_matrix[i][i];
		for(k = i; k<3; k++)
			temp_matrix[i][k] /= multiple;
		for(k = 0; k<3; k++)
			result[i][k] /= multiple;

		//for each element up from the ith diagonal element, try to make it zero.
		for(j = i - 1; j >= 0; j--)
		{
			multiple = temp_matrix[j][i];
		
			for(k = i; k < 3; k++)
				temp_matrix[j][k] -= temp_matrix[i][k] * multiple;
			for(k = 0; k < 3; k++)
				result[j][k] -= result[i][k] * multiple;
		}
	}

	return result;
}

//This method calculates the inverse of a matrix by using the gaussian elimination method with partial pivoting
double ** partialPivotInversion( double ** matrix ){
	int i = 0;
	int j = 0;
	int k = 0;
	//insert your code here.//

	double temp_matrix[3][3];

	//assign the matrix values that from the functional parameter matrix
	for(i = 0; i < 3; i ++)
		for(j = 0; j < 3; j++)
			temp_matrix[i][j] = matrix[i][j];

	//declare and allocate the result
	//initialize the result matrix as an identity matrix
	double ** result = new double*[3];
	for(i = 0; i<3; i++){
		result[i] = new double[3];
		for(j = 0; j<3; j++)
			result[i][j] = (i == j) ? 1.0 : 0.0;
	}
	
	//calculate the inverse of the matrix
	//here we do row transformation for the larger matrix of (temp_matrix, result)
	//when temp_matrix becomes identity matrix, the result matrix becomes the inverse of the original
	//temp_matrix
        
	for(i = 0; i<3; i++)
	{
		//the key feature of partial pivoting
		//find the line containing the maximun absolute value as the pivot.
		double max_fabs = 0.0;
		int max_pos = i;
		for(j = i; j<3; j++)
		{
			if(max_fabs < fabs(temp_matrix[j][i]))
			{
				max_fabs = fabs(temp_matrix[j][i]);
				max_pos = j;
			}
		}

		//swap the two lines so that the pivot position has the maximum absolute value.
		if(max_pos != i)
		{
			for(k = i; k<3; k++)
				swap(temp_matrix[max_pos][k], temp_matrix[i][k]);
			for(k = 0; k<3; k++)
				swap(result[max_pos][k], result[i][k]);
		}
		//if line i contains a zero in the diagonal place, i.e. place [i][i]
		//then we search each following lines to find an alternative line which does not contain
		//a zero in the place [alternative][i]
		int alternative = i;
		while(alternative < 3 && temp_matrix[alternative][i] == 0)
			++alternative;

		//line i and all lines after i contain zero in the column i, therefore, 
		//it can be proved that such a matrix do not have an inverse.
		if(alternative == 3)
		{
			printf("This matrix is singular and do not have an inverse!\n");
			exit(0);
		}

		//if the program runs here, it means that a non-zero alternative exists
		//if this alternative line different from line i
		//we swap each element in both rows
		if(alternative != i)
		{
			for(k = i; k<3; k++)
				swap(temp_matrix[alternative][k], temp_matrix[i][k]);
			for(k = 0; k<3; k++)
				swap(result[alternative][k], result[i][k]);
		}

		//for each row down from the ith row
                for(j = i + 1; j<3; j++)
                {
			if(temp_matrix[j][i] != 0)
			{
				//here we gurantee that matrix[i][i] is not zero
				double multiple = temp_matrix[j][i] / temp_matrix[i][i];
				for(k = i; k<3; k++)
                                	 temp_matrix[j][k] -= temp_matrix[i][k] * multiple;
				for(k = 0; k<3; k++)
					 result[j][k] -= result[i][k] * multiple;
			}
			else
			{
				//the target row has been zero, therefore, 
				//we just skip it and move on to the next row
				continue;
			}
		}
	}

	//we do the substitution for each line from the largest line to the smallest. i.e., 0
	for(i = 2; i >= 0; i--)
	{
		//try to make the ith diagonal place 1
		double multiple = temp_matrix[i][i];
		for(k = i; k<3; k++)
			temp_matrix[i][k] /= multiple;
		for(k = 0; k<3; k++)
			result[i][k] /= multiple;

		//for each element up from the ith diagonal element, try to make it zero.
		for(j = i - 1; j >= 0; j--)
		{
			multiple = temp_matrix[j][i];
		
			for(k = i; k < 3; k++)
				temp_matrix[j][k] -= temp_matrix[i][k] * multiple;
			for(k = 0; k < 3; k++)
				result[j][k] -= result[i][k] * multiple;
		}
	}

	return result;
}






/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// ###     ###     ###    ####### ###     ###    ###     ### ####### ####### ###   ###   ######   ######  
// ####   ####    #####   ####### ####    ###    ####   #### ####### ####### ###   ###  ########  ####### 
// ##### #####    #####     ###   #####   ###    ##### ##### ###       ###   ###   ### ########## ########
// ###########   ### ###    ###   ######  ###    ########### #######   ###   ######### ####  #### ###  ###
// ### ### ###   ### ###    ###   ### ### ###    ### ### ### #######   ###   ######### ###    ### ###  ###
// ###  #  ###  ###   ###   ###   ###  ######    ###  #  ### ###       ###   ######### ####  #### ###  ###
// ###     ###  #########   ###   ###   #####    ###     ### ###       ###   ###   ### ########## ########
// ###     ###  ######### ####### ###    ####    ###     ### #######   ###   ###   ###  ########  ######## 
// ###     ###  ###   ### ####### ###     ###    ###     ### #######   ###   ###   ###   ######   ######  

//execution format : [executable name] [argument]
int main(int argc, char* argv[])
{
	///usage handling
	if( argc == 1 ){
		defaultOutput();
	}

	///file format handling
	if( argc == 2 && ( strcmp(argv[1], "-format")==0 ) ){
		fileFormat();
	}

	///invert
	if( argc == 2 && ( strcmp(argv[1], "-format")!=0 ) ){
		double ** matrix = parseMatrixFile( argv[1] );
		
		//compute the simple inverse and print it
		double ** simpleInverse = simpleMatrixInversion( matrix );
		printMatrix(simpleInverse);
		
		//compute the inverse with partial pivoting and print it
		double ** partialInverse = partialPivotInversion( matrix );
		printMatrix(partialInverse);
		
		//clean up
		for(int i = 0; i<3; i++){
			delete[](simpleInverse[i]);
			delete[](partialInverse[i]);
			delete[](matrix[i]);
		}
		delete[](simpleInverse);
		delete[](partialInverse);
		delete[](matrix);
		
	}

	return 0;
}



////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////


